Contents

library(tidyverse)
library(magrittr)
library(pheatmap)
library(RColorBrewer)

1 Load data

1.1 Variance-stabilized transformed RNAseq counts

# mean vst values for each gene in each age group
counts <- read.delim(paste0(dir, "Data/vst_counts.csv"), sep = ",") %>%
  subset(select = c(ensembl_id, gene, age_group, counts)) %>%
  group_by(ensembl_id, gene, age_group) %>%
  summarize(counts = mean(counts)) %>%
  ungroup() 
## `summarise()` has grouped output by 'ensembl_id', 'gene'. You can override
## using the `.groups` argument.
head(counts)
## # A tibble: 6 × 4
##   ensembl_id      gene   age_group counts
##   <chr>           <chr>  <chr>      <dbl>
## 1 ENSG00000000003 TSPAN6 1-15        8.74
## 2 ENSG00000000003 TSPAN6 16-26       9.08
## 3 ENSG00000000003 TSPAN6 27-60       9.12
## 4 ENSG00000000003 TSPAN6 61-85       9.01
## 5 ENSG00000000003 TSPAN6 86-96       9.46
## 6 ENSG00000000419 DPM1   1-15        9.48

1.2 DE genes

# DE genes
#age <- data.frame('transition' = c('fc_16-26_27-60', 'fc_61-85_86-96'), 
#                  'DE' = c('young', 'old'))
DE_genes <- read.delim(paste0(dir, "Data/DE_var_p_n_200.csv"), sep = ",") %>%
  #inner_join(age) %>%
  subset(select = c(gene, transition))
head(DE_genes)
##    gene     transition
## 1   A2M  fc_1-15_16-26
## 2   A2M fc_16-26_27-60
## 3  AARD  fc_1-15_16-26
## 4 ABCA3 fc_27-60_61-85
## 5 ABCA5 fc_16-26_27-60
## 6 ABCA5 fc_61-85_86-96

2 Expression Heatmap

2.1 Group 1 vs. Group 2

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  filter(transition == 'fc_1-15_16-26') %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id') %>%
  subset(select = -transition)
## Joining, by = "gene"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

2.2 Group 2 vs. Group 3

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  filter(transition == 'fc_16-26_27-60') %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id') %>%
  subset(select = -transition)
## Joining, by = "gene"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

out <- pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

gene_order <- rownames(z_scores[out$tree_col[["order"]], ])
gene_order <- c(gene_order[27:162], gene_order[1:26])
pheatmap(t(z_scores[gene_order, ]), cluster_rows = FALSE, cluster_cols = FALSE, show_colnames = FALSE, fontsize = 22)

2.3 Group 3 vs. Group 4

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  filter(transition == 'fc_27-60_61-85') %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id') %>%
  subset(select = -transition)
## Joining, by = "gene"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

2.4 Group 4 vs. Group 5

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  filter(transition == 'fc_61-85_86-96') %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id') %>%
  subset(select = -transition)
## Joining, by = "gene"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

out <- pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

gene_order <- rownames(z_scores[out$tree_col[["order"]], ])
pheatmap(t(z_scores[gene_order, ]), cluster_rows = FALSE, cluster_cols = FALSE, show_colnames = TRUE, 
         fontsize = 22, fontsize_col = 12)

2.5 Group1-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group1.csv"), sep = ",") %>%
  subset(select = c(ensembl_id))

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "ensembl_id"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

2.6 Group4-5-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group4_5.csv"), sep = ",") %>%
  subset(select = c(gene))

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "gene"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

2.7 Group5-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group5.csv"), sep = ",") %>%
  subset(select = c(ensembl_id))

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "ensembl_id"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

2.8 Upregulated Group5-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group5_up.csv"), sep = ",") %>%
  subset(select = c(ensembl_id))

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "ensembl_id"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")
breakslist = seq(-1.5, 1.5, by = 0.1)

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 30,
         breaks = breakslist, color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(breakslist)))

2.9 Downregulated Group5-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group5_down.csv"), sep = ",") %>%
  subset(select = c(ensembl_id))

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "ensembl_id"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

breakslist = seq(-1.5, 1.5, by = 0.1)

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 30, 
         breaks = breakslist, color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(breakslist)))

2.10 Upregulated Group45-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group45_up.csv"), sep = ",") %>%
  subset(select = c(gene))

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "gene"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

2.11 Downregulated Group45-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group45_down.csv"), sep = ",") %>%
  subset(select = c(gene))

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "gene"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

2.12 Upregulated Group4-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group4_up.csv"), sep = ",") %>%
  subset(select = c(gene))

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "gene"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

2.13 Downregulated Group4-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group4_down.csv"), sep = ",") %>%
  subset(select = c(gene))

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "gene"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 22)

2.14 Upregulated Group1-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group1_up.csv"), sep = ",") %>%
  subset(select = c(ensembl_id))

z_scores <- inner_join(counts, DE_genes) %>%
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "ensembl_id"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")
breakslist = seq(-1.5, 1.5, by = 0.1)

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 30, 
         breaks = breakslist, color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(breakslist)))

2.15 Downregulated Group1-specific genes

DE_genes <- read.delim(paste0(dir, "Data/DE_Group1_down.csv"), sep = ",") %>%
  subset(select = c(ensembl_id))

z_scores <- inner_join(counts, DE_genes) %>% 
  subset(select = -c(gene)) %>%
  group_by(ensembl_id) %>%
  mutate(counts = scale(counts)) %>%
  spread(key = age_group, value = counts) %>%
  column_to_rownames(var = 'ensembl_id')
## Joining, by = "ensembl_id"
colnames(z_scores) <- c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5")
breakslist = seq(-1.5, 1.5, by = 0.1)

pheatmap(t(z_scores), cluster_rows = FALSE, show_colnames = FALSE, fontsize = 30, 
         breaks = breakslist, color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(breakslist)))